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Abstract 

Background: The rice white tip nematode Aphelenchoides besseyi, a devastating nematode whose genome has not been 
sequenced, is distributed widely throughout almost all the rice-growing regions of the world. The aims of the present study 
were to define the transcriptome of A. besseyi and to identify parasite-related, mortality-related or host resistance- 
overcoming genes in this nematode. 

Methodology and Principal Findings: Using Solexa/lllumina sequencing, we profiled the transcriptome of mixed-stage 
populations of ;A. besseyi. A total of 51,270 transcripts without gaps were produced based on high-quality clean reads. Of all 
the A. besseyi transcripts, 9,132 KEGG Orthology assignments were annotated. Carbohydrate-active enzymes of glycoside 
hydrolases (GHs), glycosyltransferases (GTs), carbohydrate esterases (CEs) and carbohydrate-binding modules (CBMs) were 
identified. The presence of the A. besseyi GH45 cellulase gene was verified by in situ hybridization. Given that 13 unique A. 
besseyi potential effector genes were identified from 41 candidate effector homologs, further studies of these homologs are 
merited. Finally, comparative analyses were conducted between A. besseyi contigs and Caenorhabditis elegans genes to look 
for orthologs of RNAi phenotypes, neuropeptides and peptidases. 

Conclusions and Significance: The present results provide comprehensive insight into the genetic makeup of A. besseyi. 
Many of this species' genes are parasite related, nematode mortality-related or necessary to overcome host resistance. The 
generated transcriptome dataset of A. besseyi reported here lays the foundation for further studies of the molecular 
mechanisms related to parasitism and facilitates the development of new control strategies for this species. 
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Introduction 

Parasitic nematodes that prey on plants cause devastating losses 
in agricultural production and contribute to a significant reduction 
in yield. It has been estimated that plant parasitic nematodes cause 
annual losses in crop production valued at US$125 billion [1]. The 
rice white tip nematode Aphelenchoides besseyi Christie, 1942, which 
is widely distributed throughout almost all the rice growing regions 
of the world, is considered to be a major contributor to the seed- 
borne pathogens of rice [2^] . A. besseyi is a global burden, causing 
losses of up to 50% in upland rice in Brazil [5] and 10-50% in 
flooded rice in China [6]. 

As a seed-borne nematode, A. besseyi can survive on stored grain 
in anhydrobiosis for several years [7] . This species is a facultative 
migratory ecto- and endoparasite of the leaf or young tissue of rice, 
causing whitening of the top of the leaf, which later becomes 
necrotic, and a crinkling and distortion of the flag leaf enclosing 
the panicle, which then dies off and disintegrates [4,8]. Moreover, 
A. besseyi is polyphagous and damages a wide range of host plants. 



including many important crop species, such as strawberry 
(Fragaria grandiflord), sweet corn [^ea mays), sweet potato [Dioscorea 
escuknta), bird's nest fern {Asplenium nidus) and many vegetables and 
grasses [9-1 1]. 

Based on its phylogenetic tree [12], A. besseyi belongs to clade 
lOB, whereas most other major plant-parasitic nematodes belong 
to clade 12. Only a high-quality draft genome sequence from 
Bursaphelenchus xybphilus is available for the investigation of the 
biological basis of clade 10 nematodes [13]. The investigation of A 
xylophilus genes linked to key biological processes and parasitism 
will provide guidance for further studies of A. besseyi and the other 
clade 10 nematodes. Although the model free-living nematodes 
Caenorhabditis spp. (clade 9A), the human parasite Brugia malayi 
(clade 8B) and several plant-parasitic nematodes (clade 12) have 
been subjected to molecular analysis [14], very little is known 
about the genomic background of ^. besseyi. However, recently, a 
fatty acid and retinoid-binding protein gene [Ab-far-l) was cloned 
from A. besseyi based on EST sequencing [15]. In addition, an 
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exciting small-scale transcriptome project for A. besseyi was recently 
published [16,17]. 

Next-generation sequencing technology, such as the Solexa/ 
lUumina, 454 (Roche) and ABI-SOLiD platforms, has dramati- 
cally improved the efficiency of gene discovery [18]. Using next- 
generation sequencing technology, several nematode transcrip- 
tomes have been sequenced, including Trichuris suis, Mecator 
americanus, B. malayi, Trichostrongyltis colubriformis, Praplenchus cojfeae 
and P. thomei [1 4, 1 9-24] . Whereas genomic analyses of A xykphilus 
reveal the genetic potential of clade 1 0 nematodes, transcriptomic 
studies of nematodes will enable the regulator)' networks 
controlling the expression of parasitism-related genes in space 
and at dilferent stages to be deciphered. 

In the present study, Solexa/IUumina sequencing technology 
was used to explore and functionally annotate the transcriptome of 
a mixed-stage population of ^. hes.seyi. The results provide insights 
into the genetic background of ^. besseyi. The transcriptomic data 
presented here serve as a valuable resource that will help deepen 
our understanding of the parasitism of this nematode. Under- 
standing the basic biology of plant-parasitic nematodes, especially 
the parasitic mechanisms at the molecular level, is essential for the 
development of effective control strategies. 

Materials and Methods 

Nematode Material 

The nematode A. besseyi (NCBI BioSample accession No. - 
SAMN02420038), derived from a single female isolated from an 
infected rice grain in eastern China, was maintained by serial 

subculture on a Bottytis cinerea mat grown on sterilized barley [25] 
and kept in the dark at 25°C. Nematodes freshly extracted from 
sterilized barley were used for RNA extraction. Eggs, juvenile 
stages (J2^4) and adults were mixed together at a ratio of 1:2:3. 

RNA Extraction, cDNA Library Construction and lllumina 
Sequencing 

Total RNA was extracted from mixed-stage nematodes using 
TRIzol regent (Invitrogen, Carlsbad, CA, USA) [26] and then 
cleaned using an RNeasy Minikit column (Qiagen, Valencia, CA, 
USA) according to the manufacturers' instructions. The extracted 
RNA was assessed for quality and quantified using a BioPhotom- 
eter D3() (Eppendorf, Hamburg, Germany). 

Two micrograms of total RNA was used to extract polyadenyl- 
ated RNA and construct strand-specific RNAseq libraries accord- 
ing to the lllumina protocol (directional mRNASeq sample 
preparation part #15018460 Rev. A, October 2010; lllumina, 
San Diego, CA, USA) (detailed in Text SI in File SI). The library 
quality was validated with a Bioanalyzer 2100 (Agilent, CA, USA). 
Each sample was sequenced for 80 cycles on one lane of the 
lUumina GAIIX platform (lllumina, San Diego, CA, USA) with a 
2x100 bp module by LC Sciences (Houston, TX, USA). The 
dataset generated by deep-sequencing platforms was deposited in 
the NIH short-read archive (accession No. SRR1040470). 

Data Filtering and De Novo Assembly 

Low-quality reads with ambiguous sequences "N" were 
removed according to a sliding \\in(l()\\" mc'tliod using a threshold 
of 35 bp. To examine the coverage of the sequences, 500,000 
randomly selected reads were compared with the NT (non- 
redundant nucleotide sequences in the NCBI) database using 
BLAST, and a typical cutoff expectation value (E-value) was less 
than le-10. The trinity (trinityrnaseq_r20 13-02-25) [27] and 
paired-end methods were used to generate non-redundant 
unigenes via the de novo assembly of the clean reads. The quality 



of the assembly was critically assessed by LC Sciences (Houston, 
TX, USA) before subsequent analysis. 

The expression level of each transcript was measured as the 

number of clean reads mapped to its sequence. The mapped clean 
read number was normalized to RPKM (reads per kilobase of 
exon model per million mapped reads) with Bowtie 0.12.8 (single- 
end mapping) and adjusted with a normalized factor [28]. 

To determine how many unigenes in total were available in our 
transcriptome and another recently published A. besseyi transcrip- 
tome [16], we compared the unigenes with the 5,804 MIRA 
assembly sequences downloaded from nematode.net [29] using 
BLASTN widi an E-value cutoff of le-6. 

Annotation and KEGG Pathway Determination 

Unigene sequence directions were determined by searches 
against the NR (non-redundant protein sequences in the NCBI), 
Swiss-Prot, TrEMBL, Cdd, Pfam and KOG databases (E- 
value<le-5). 

KEGG mapping was used to determine the metabolic pathways 
of metabolism, genetic information processing, environmental 
information processing, cellular processes, organismal systems and 
etc. [30]. Sequences were submitted to the KEGG Automatic 
Annotation Server (KAAS) [31] to enrich the pathway annotation. 

Annotation of Carbohydrate-Active Enzymes 

The CAZymes Analysis Toolkit [32] was used to detect A. besseyi 
carbohydrate-active enzymes (CAZymes) according to a previous 
report [13]. Expansin-like proteins, pectinases and arabinanases 
were detected with a BLAST search using con; modules of kn(nvn 
protein sequences as queries. BLASTP was used to predict protein 
functions by searching against known protein modules and 
catalytic sites downloaded from NCBI's Conserved Domain 
Database service. Multiple sequence alignments were performed 
using Clustal W, and phylogenetic trees were constructed based on 
those sequences using the MEGA 5.05 program. 

Effector Candidate Identification 

Known putative effectors and selected annotated secreted 
protein sequences from different plant-parasitic nematodes 
(Aphetenchus avenae, B. mucromtus, B. xylophilus, Ditylenchus afiicams, 
Ghbodera mexicana, G. pallida, G. rostochiensis, G. tabacum, Heterodera 
glycines, H. schachtii, M. arenaria, M. chitwoodi, M. hapla, M. incognita, 
M. javanica, P. cojfeae, P. penetrans, Rotylenchulus reniformis, Radopholus 
similis and Xiphinema index) [33] were downloaded from GenBank. 
These sequences were searched against the A. besseyi contigs with 
TBLASTN to identify potential homologs [34] . An E-value cutoff 
of le-5 was used to identify significant matches. SignalP 4.1 [35] 
and TMHMM 2.0 [36] were used to identify secreted proteins as 
previously described [13]. 

RNAi Phenotype iVlatches to C. elegans Genes 

Seventy-five C. elegans in vitro RNA interference (RNAi)-related 
genes were identified and retrieved from WormBase (version 

WS236). The A. besseyi contigs were compared with C. elegans genes 
to identify orthologs of RNAi phenotypes [37]. 

Neuropeptides and Peptidase Detection 

One hundred thirteen C. elegans neuropeptide genes were 
identified and retrieved from WormBase (\ ersion WS236). The A. 
besseyi contigs were compared with C. elegans genes to identify 
orthologs with high identity to neuropeptides by BLAST. 
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■ Ascaris suum 

■ Loa loa 

Caenorhabditis elegans 

■ Brugia malayi 

■ Caenorhabditis brenneri 

■ Caenorhabditis remanei 

■ Caenorhabditis briggsae 

■ Bursapheienchus xylophiius 
I Haemonchus contortus 

■ Heterodera glycines 

Figure 1. The ten most common NR hierarchies identified in the A. besseyi sequences. 
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The MEROPS database 9.8 was used to identify putative A. 
besseyi peptidases [38] . The candidates were derived and examined 
for similarity (E-vaIue< Ie-10) to MEROPS proteins. 

In Situ Hybridization of GH45 Cellulases 

In situ hybridization was performed on A. besseyi essentially as 
described by De Boer et al. [39]. A digoxigenin-labeled A. besseyi 
GH45 cellulase DNA probe was synthesized using the cDNA 
sequence following the general protocol of the Roche DIG DNA 
Labeling Kit (Roche, Mannheim, Germany, cat. 
no. 11175033910). The specificity of the probe was verified by 
Southern blotting. The hybridization signals were detected by 
incubation with the DIG High Prime DNA Labeling and 
Detection Starter Kit I (Roche, Mannheim, Germany, cat. 
no. 11745832910) and photographed under a Bx51 microscope 
(Olympus, Tokyo, Japan). 

Results 

lllumina Sequencing and De Novo Assembly 

A total of 46,826,350 raw reads were obtained by Solexa 
(lUumina GAIIX) RNA paired-end sequencing, and 36,905,372 
high-quality clean reads were obtained after a strict filtering 
process. The average length of the clean reads was 89.89 bp. The 
ratio of clean reads was 78.81%. No RNA contamination from B. 
cinerea or other organisms was found based on 500,000 randomly 
selected reads from the NT database. Using the clean reads. 
Trinity produced 51,270 transcripts with length >200 bp. The 
average length and N50 of these transcripts were 1,241 bp and 
1,857 bp, respectively. A total of 35,180 unigenes with an average 
length of 1,050 bp and N50 of 1,626 bp was generated after 
further clustering and assembly. Of these unigenes, 21,882 
unigenes (62.2%) were > 500 bp and 13,259 unigenes (37.7%) 
were > 1,000 bp (Figure SI in File SI). The frequency distribution 
of the GO content is shown in Figure S2 in File SI. 

Recently, another A. besseyi transcriptome including 5,804 
MIRA assembly sequences was published [16]. Compared with 
this transcriptome (detailed in the supplementary dataset in File 



S2), 77.8% (27,380 out of 35,180) of the unigenes in our 
transcriptome were unique. Only 2.2% (125 out of 5,804) of the 
MIRA assembly-yielded sequences were not found in our dataset 
(detailed in Text S2 in File SI). Because the E-value cutoff of le-6 
was used, 7,800 (22.2% of the total 35,180) of the unigenes in our 
transcriptome matched 5,679 (97.8% of the total 5,804) MIRA 
assembly sequences. Therefore, at least 33,184 unigenes are now 
available in the two datasets combined. 

The distribution of transcript values indicated that all transcripts 
were expressed, and 95.69% of the reads were mapped (Table SI 
in File SI). Of the 51,270 transcripts, die RPKM values of the 
majority of the transcripts varied from 1 to 12 (Table S2 and 
Figure S3 in File SI). The distribution of RPKM values indicated 
that most genes were expressed at low levels. 

Annotation 

Unigenes were determined by performing BLASTX searches 
against the NR protein database. A total of 17,338 orthologs were 
identified. Of the 35,180 unigenes, 88.03% were annotated as the 
most common ten organisms (Figure 1). The majority of the 
annotated A. besseyi (clade lOB) sequences corresponded to known 
nucleotide sequences of nematodes belonging to clades 2A 
[Trichinella spiralis), 8B [Ascaris suum, Loa loa and B. malayi), 9A 
[Caenorhabditis spp.) and lOB [B. xylophiius). Other protein databases, 
such as Swiss-Prot, TrEMBL, Cdd, Pfam and KOG, were used to 
annotate these unigenes (detailed in Table S3 in File SI). A search 
for Gene Ontology (GO) identifiers for these unique sequences was 
performed, and 14,064 GO terms were retrieved. The most 
abundant GO terms in the dataset are shown in Figure 2. 

Pathway Determination 

KAAS was used to annotate the EST sequences of the animal 
nematodes [B. malayi and A. suum), the free-living nematode [C. 
elegans), A. besseyi and other plant-parasitic nematode [B. xylophiius 
and M. incognita) transcripts. Of all the A. besseyi transcripts, 9,132 
KEGG Orthology (KO) assignments were annotated. Two Venn 
diagrams were drawn to show the intersections of the nematode 
KO assignments (Figure 3). The diagram of C. elegans, A. besseyi and 
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biological prioces 



cellular component 



Figure 2. The most abundant Gene Ontology terms present in tKie dataset for the cellular component, molecular function and 
biological process categories. 
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the other plant parasitic nematodes indicated 8,705 A. besseji KO 
assignment hits in one or more of the nematodes, of which 427 
KO assignments were unique (Figure 3A). Another diagram of C. 
elegans, A. besseji and the animal-parasitic nematode A. suum also 
illustrates 429 unique A. hesseji KO assignments (Figure 3B). 



To identify the pathways in which the unigenes participated, we 
mapped the A. bessefi unigenes to the canonical reference pathways 
in the KEGG database. A total of 9,133 unigenes were annotated 
based on a BLASTX search of the KEGG database; 753 enzymes 
and 306 KEGG pathways were predicted (Table 1, Tables S4— S8 
in File SI). Of these KEGG pathways, 121 pathways were 
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Figure 3. The intersections of lists of nematode KEGG Orthology (KO) assignments. A: Assignment of C. elegans, A. besseyi and other plant- 
parasitic nematodes. B: Assignment of C. elegans, A. besseyi and the animal-parasitic nematode A. suum. Ab: A. besseyi, Bx: 6. xyloptiilus, Mi: M. 
incognita, Ce: C. elegans, Bm: 6. malayi. As: A. suum. 
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consistent with the corresponding pathways in C. elegans. In 
addition to the most common 5 pathways (Table 2), other 
important patliways, including ubiquitin-mediated proteolysis (72 
members), glutathione metabolism (39 members), glycolysis/ 
gluconeogenesis (38 members), glycerophospholipid metabohsm 
(38 members), drug metabolism (35 members) and metabolism of 
xenobiotics by cytochrome (32 members), were found. 

CAZymes and Cell Wall Degradation 

Some CAZymes are known to be secreted by sedentary as well 
as migratory plant parasitic nematodes and have been predicted to 
be involved in nematode migration in the plant and feeding from 
host cells [40]. According to the Carbohydrate- Active enZYmes 
Database, CAZymes are grouped into the following classes: 
glycoside hydrolases (GHs), glycosyltransferases (GTs), polysac- 
charide lyases (PLs), carbohydrate esterases (CEs) and carbohy- 
drate-binding modules (CBMs) [24]. A summary of the CAZymes 
detected in A. besseji is provided in Table 3. In all, 524 contigs were 
assigned to the following four classes of CAZymes: 136 contigs 
with GH domains to 10 families of GHs (Table S9 in File SI), 328 
contigs to 28 families of GTs (Table SIO in File SI), 42 contigs to 2 
families of CEs (Table SI 1 in File SI) and 18 contigs to 3 families 
of CBMs (Table S12 in File SI). Apparently, PLs are not present in 
any of the contigs in this transcriptome. 

With regard to plant cell wall degradation-related CAZymes, we 
first noted that 4 contigs were identified as homologs of GH45 
ceUulases previously identified in B. xylophilus and another A. besseyi 
transcriptome [41]. Multiple sequence alignments of A. besseyi 
GH45 ceUulase with B. xylophilus homologs and those of 
Ascomycota and Basidiomycota fungi revealed that the motifs 
were almost identical in all proteins (Figure 4A). The GH45 
ceUulase maximum likelihood tree (Figure 4B) showed four major 
clusters. Because no other nematode GH45 ceUulases could be 
found, B. xylophilus and B. mucronatus were the only two plant 
parasitic nematodes used in the phylogenetic tree study in addition 
to A. besseyi. 

In situ hybridization was used to determine the tissue specificity 
of A. besseyi GH45 ceUulase gene transcription. The anti-sense 
probe stained the esophageal gland ceU area (Figure 5). The gland 
region was strongly stained by probes against the A. besseyi GH45 
ceUulase gene, whereas the sense probes gave no signal (Figure 5E). 
The A. besseyi GH45 ceUulase gene was expressed at a lower level 
in males (Figure 5A and B) than in females (Figure 5C and D). 

Another plant cell wall degradation-related CAZyme famUy 
secreted by plant parasitic nematodes is GH5 ceUulase, which has 
been weU studied in the plant-parasitic nematodes of the order 
Tylenchida (clade 12, including most of the major plant-sedentary 
endoparasitic nematodes such as Meloidogyne spp., Heterodera spp. 



Table 1. The classification A. besseyi pathways. 







Pathways 


Members 


Metabolism 


121 


2077 


Genetic information Processing 


22 


868 


Environmental Information Processing 


25 


478 


Cellular Processes 


17 


583 


Organismal Systems 


60 


957 


Human Diseases 


61 


1396 


Total 


306 


4489 
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Table 2. The most common 5 A. besseyi pathways consistent 
with C. elegans. 



Pathway Pathway type Unigene num 

celOl 100 Metabolic pathways 602 

cel03010 Ribosome 119 

cel04141 Protein processing in endoplasmic reticulum 118 

cel00190 Oxidative phosphorylation 106 

cel0301 3 RNA transport 99 

doi:1 0.1 371/journal.pone.0091 591 .t002 



Table 3. Carbohydrate-active enzymes identified in A. besseyi. 





CAZy enzyme classes 


CAZy family 


A. besseyi con\\gi 


Glycoside hydrolases 


10 


136 


Glycosyl transferases 


28 


328 


Carbohydrate esterases 


2 


42 


Polysaccharide lyases 


0 


0 


Carbohydrate-binding modules 


3 


18 


Total 


43 


524 



doi:1 0.1 371 /journal.pone.0091 591 .t003 



and Globodera spp. and the migratory parasitic species P. thornei and 
P. cqffeae) [14,24]. However, GH5 ceUulases were apparenfly 
absent from the A. besseyi transcriptome. Furthermore, GH28 
polygalacturonase, GH30 xylanase, GH43 arabinase, GH53 
arabinogalactan galactosidase, PL3 pectate lyase and expansin- 
like protein were absent. 

In addition to plant cell waU degradation-related CAZymes, we 
also noted that A. besseyi possessed GH16 and GH20 CAZymes, 
which may be involved in fungal cell wall degradation (Table S9 in 
FUe SI). GH16, which degrades beta-l,3-glucan, a core compo- 
nent of the fungal cell waU, was present in this A. besseyi 
transcriptome in remarkable abundance (a total of 90 contigs). 
In contrast, GH20, which is involved in chitin degradation, was 
detected in low abundance (2 contigs). 

Candidate Effector Identification 

Numerous effectors with various functions in the plant cell are 
secreted into the plant by nematodes [33]. Unlike sedentary or 
migratory endoparasites, A. besseyi is a migratory ecto- and 
endoparasite nematode that does not rely on biotrophy. Candidate 
A. besseyi parasitism genes and effector homologs were identified to 
estabhsh the relationships between the parasite strategy and 
parasitism genes. AU known and annotated putatively secreted 
proteins from plant-parasitic nematodes were compared with the 
A. besseji contigs to identify possible homologs. 

In addition to the CAZymes described above, 41 candidate 
effector homologs were identified (Table S13 in File SI). Thirteen 
unique A. besseyi potential effector genes were identified (Table 4) 
from the 41 candidate effector homologs based on the defining 
characteristic that effector genes that encode proteins have a 
putative signal peptide but lack a transmembrane helix and 
have an open reading frame of at least thirty amino acids after 
the predicted signal peptide cleavage site [42-44]. In contrast 
to sedentary endoparasites, no plant-hormone mimics or 
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Figure 4. Alignment and polygenetic tree of the A. besseyi GH45 cellulase amino acid sequence with its homologs identified from 
the NCBI database. A: The alignment shows the comparison between A. besseyi GH45 cellulase and its homologs. Active site residues are marked 
with red asterisks. The GH45 cellulase consensus sequences are boxed with a red line and highlighted with a magnifying lens. Asterisks, double dots 
and single dots denote fully conserved, strongly conserved and weakly conserved amino acid residues, respectively. B: The maximum likelihood tree 
shows the relationship between A besseyi GH45 cellulase and homologs. One thousand bootstrap replicates were performed, and the node labels 
represent the percentage of bootstrap support. 
doi:1 0.1 371/journal.pone.0091 591 .g004 



nodulation-related effector homologs were found in the A. 
transcriptome. 

Because A. besseyi does not establisli a long-term feeding 
relationship with the living cells of their hosts, it was necessary 
to confirm the 13 potential effectors by functional classification. 
Similar to B. xylophilus [13], the majority (7 out of the 13) of^. 



besseyi candidate effectors were also present in other migratory 
plant nematodes, C. elegans and even animal-parasitic nematodes. 
Cathepsin-like proteinases were also found in another migratory 
nematode, R. similis [45]. Chitinase, an egg-hatching-related 
protein, was found in Caenorhabditis spp. and L. loa. Glutathione 
peroxidase was found in Caenorhabditis spp., A. smm and Haemonchus 
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genes are previously assumed to be involved in the process of 
syncytia or giant cell formation, these findings suggest that the 
roles of these genes may be different in A. besseyi. 

The Presence of Genes Involved in the RNAi Pathway 

RNAi is an effective method for gene function identification and 
will facilitate the identification of potential targets for developing 
potent antinematodal drugs against A. hesseji. RNAi has been 
successfully applied in both sedentary and migratory plant- 
parasitic nematodes, including more than twenty plant-parasitic 
nematodes [46-50] . Recendy, RNAi was successfully applied to A. 
besseyi [15]. 

A comparative analysis was conducted with the A. besseyi contigs 
and C. elegans genes to look for RNAi phenotype orthologs. The 
2,899 A. besseyi contigs homologous to C. elegans matched 57 RNAi 
pathway effectors. The RNAi-related proteins were separated into 
the following five core functional groups: small RNA biosynthesis, 
dsRNA uptake and spreading, argonautes (AGOs) and RNA- 
induced silencing complex (RISC), RNAi inhibitors and nuclear 
effectors (Table 5 and Tables S14-S18 in File SI). Apparently, A. 
besseyi encodes more predicted RNAi pathway effectors than B. 
xjhphilus [13] and M. incognita [36]. 

Neuropeptide Detection 

Neuropeptides play important roles in the modulation of 
synaptic activity and the development of the nervous system and 
may also function as primary neurotransmitters [24,51]. Accord- 
ing to WormBook, 113 C. elegans neuropeptide genes, which 
encode over 250 distinct neuropeptides, have been identified [52]. 

With the sequencing of the A. besseyi transcriptome, 49 
neuropeptide genes were identified (Table 6). Of these, 8 genes 
encode insulin-like peptides (Table S19 in File SI), 17 genes 
encode FMRFamide-related peptides (Table S20 in File SI) and 
24 genes encode non-insulin or noii-FLP peptides (Table S21 in 
File SI). 

Peptidase Detection 

The MEROPS database classifies peptidases into seven catalytic 
types; aspartic, glutamic, cysteine, isoleucine, metallo, serine and 
threonine peptidases. In our analysis, 73 peptidase genes (Table 7 



Table 4. Potential effectors identified in A. besseyi. 





Effector 


Function 


Contigs with hit 


Bit score 


Best E-value 


4D01 


unl<nown 


21 


135 


2e-32 


5G05 


unl<nown 


6 


175 


5e-44 


6F06 


unl<nown 


56 


66 


5e-ll 


Acid phosphatase 


digestion 


30 


324 


le-88 


Calreticulin 


calcium signaling 


24 


617 


e-177 


Cathepsin 


protein degradation 


2 


449 


e -126 


Chitinase 


egg hatcfiing 


7 


125 


9e-29 


C-type lectin 


unl<nown 


3 


48 


4e-05 


FAR 


binding of host fatty acids; reducing defense response 


12 


177 


4e-45 


Glutathione peroxidase 


detoxification of ROS 


10 


360 


e-100 


PDI 


protein folding 


4 


706 


0 


Putative esophageal gland cell secretory 
protein 3 


unknown 


10 


135 


2e-32 


SXP-RAL2 


unl<nown 


2 


105 


le-23 
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Figure 5. Localization of A. besseyi GHK cellulase mRNA by in 
situ hybridization using digoxigenin-labeled antisense or sense 
probes. A and B: Antisense probe in a male nematode. C and D: 
Antisense probe in a female nematode. E: Control sense probe in a 
female nematode. The scale bars represent 20 jim. G: esophageal 
glands; S: stylet; M: metacarpus. 
doi:1 0.1 371 /journal.pone.0091 591 .g005 



contortus. FAR, which was also found in an A. besseyi EST project, 
was found in Caenorhabditis spp. and H. contortus. PDI was found in 
Caenorhabditis spp., B. malayi and L. loa. The other 6 effector 
homologs had unknown functions. Given that all the identified 
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Table 5. Presence of genes involved in the RNAi pathway. 







Number of RNAI effector proteins 






>1, /7e55ey/ contigs* 




A. besseyi 


C. elegans 


B. xylophilus 


M. incognita 






clade 10B 


clade 9A 


clade 10B 


clade 12B 




Small RNA biosynthesis 


7 


9 


8 


7 


112 


dsRNA uptal<e and spreading 


6 


11 


6 


4 


84 


AGOs and RISC 


27 


30 


12 


9 


2212 


RNAi inhibitors 


5 


10 


2 


2 


73 


Nuclear effectors 


12 


15 


9 


6 


418 


Total 


57 


75 


37 


28 


2899 



*(bitscore>40 and E-value<e-4). 
doi:l 0.1 371 /journal.pone.0091 591 .t005 



and Table S22 in File SI) were identified in A. besseyi. Proteins 
belonging to the COIA (41 contigs), 102 (89 contigs), M12A (43 
contigs) and S09X (110 contigs) families were particularly 
expanded in A. besseyi. Peptidase families involved in extracellular 
digestion and lysosomal activities were also particularly expanded. 

Discussion 

This study provides the first detailed description of the 
transcriptome of A. besseyi, a rice-parasitic nematode without a 
sequenced genome. We detected and annotated the repertoire of 
candidate CAZymes in the A. besseyi transcriptome. Our dataset 
extends the knowledge of the arsenal of enzymes in A. besseyi. 
CeUulases are classified into 14 glycoside hydrolase families — 
GH5, 6, 7, 8, 9, 10, 12, 26, 44, 45, 48, 51, 61 and 74 [53] —most 
of which are typically produced by plant pathogens and plant 
parasites [54]. GH45 ceUulases are present in Bursaphelenchus 
species and A. besseyi [13,16]. The discovery of these genes led to 
the hypothesis that they were acquired from fungi via horizontal 
gene transfer (HGT) [40]. 

Although a GH5 ceUulase was found in A. fmgmiae [55], no 
transcripts similar to GH5 ceUulases were present in the A. besseyi 
transcriptome. To confirm this result, a TBLASTN search of ^. 
fragariae GH5 ceUulase protein sequences (NCBI accession 
Nos. AFD33557, AFD33558 and AFI63769) against the A. besseyi 
contigs was performed, but no hits were found with an E-value 
cutolf of 1 e-4. It is worth noting that a small-scale transcriptome 
project for A. besseyi has been published [16]. GH45 ceUulase was 
found in this transcriptome. However, no transcripts similar to 
GH5 ceUulases were identified in the smaU-scale transcriptome. 
These results are consistent with our study. 

In our work, a remarkable abundance of GH16 (endo-1, 3-beta- 
glucanase), which is thought to be important for the degradation of 



Table 6. Neuropeptide genes in A. besseyi. 





Neuropeptide Genes 


Contigs with hit E-value level 


Insulin-like peptides 8 


25 


0.003-8e-015 


FMRFamide-related peptides 17 


54 


5e-005-6e-028 


Non-insulln, non-FLP 24 
peptides 


288 


3e-005-2e-032 


Total 49 


367 





Table 7. Summary of peptidases 


In A. besseyi. 




Catalytic types 


Family number 


Best E value 


Aspartic 


2 


1.50E-106 


Cysteine 


19 


1.60E-108 


Isoleucine 


6 


l.lOE-55 


Metallo 


26 


1.50E-148 


Serine 


16 


9.60E-100 


Threonine 


3 


5.50E-202 


Unassigned 


1 


2.10E-46 


Total 


73 




doi:1 0.1 371 /journal.pone.0a91 591 .t007 
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fungal ceU walls, was detected. This enzyme has been identified in 
various fungivorous and plant-parasitic nematodes, first in 
Bursaphelenchus species [56] and more recently in A. avenae [57]. 
Because it is a facultative fungal feeder, it makes sense that A. 
besseyi harbors considerable GH16. This result supports the 
possibUity that these fungi-consuming nematodes acquired endo- 
1,3-fi-glucanase genes from bacteria to obtain their fungal feeding 
abUity and subsequently acquired ceUulase genes from fungi, 
which permitted them to parasitize plants [56,58]. 

Of particular interest was the detection of candidate effector 
homologs such as cytokinins, ubiquitin extension protein, 
calreticulin, 10A06 and NULGla [59] (Table S13 in FUe SI), 
although these genes were thought to be involved in the process of 
either syncytia or giant cell formation [47,60-63]. Because A. 
besseyi does not generate specific feeding cells, the roles of these 
genes may indicate a dilferent function. In addition, the effector 
identification of the small-scale transcriptome project for ^. besseyi 
indicated that almost aU (12 out of 13) of the candidate sequences 
had matches in Caenorhabditis elegans and in other nematodes [16]. 
These results suggested that whether the identified A. besseyi 
homologs reaUy play the same roles as the effectors of root knot 
and cyst nematodes in parasitism is an open question. Obviously, 
further studies of these homologs are merited. 

The present results provide insights into the genetic background 
of A. besseyi. These insights will be useful in future investigations of 
nematode parasitism the development of techniques for nematode 
control and plant protection. 



PLOS ONE I www.plosone.org 



8 



March 2014 | Volume 9 | Issue 3 | e91591 



White Tip Nematode Transcriptome 



Supporting Information 

File SI Supporting Figures, Tables, and Text. Figure SI, 
Length distribution of All-Unigene. Figure S2, GC content 
frequencies distribution. Figure S3, RPKM distribution of 

transcript. Table SI, Expressed transcripts. Table S2, RPKM of 
the A. besseyi transcripts. Table S3, The annotation of A. besseyi 
unigenes. Table S4, Metabolism related pathways. Table S5, 
Genetic information processing related pathways. Table S6, 
Environmental information processing related pathways. Table 
S7, Cellular processes related pathways. Table S8, Organismal 
systems related pathways. Table S9, Glycoside hydrolases 
identified in A. hesseyi. Table SIO, Glycosyl transferases identified 
in A. hesseyi. Table Sll, Carbohydrate esterases identified in A. 
besseyi. Table SI 2, Carbohydrate binding modules identified in ^. 
besseyi. Table SI 3, Candidate effectors or potential parasitism 
genes identified in A. besseyi. Table SI 4, Small RNA biosynthesis 
proteins. Table SI 5, dsRNA uptake and spreading effectors. Table 
SI 6, Argonautes and RNA-induced silencing complex compo- 

References 

1 . Chitwood DJ (2003) Research on plant-parasitic nematode biology conducted by 
the United States Department of .^^iculture-Agricultural Research Service. Pest 
Manage Sci 59; 748-753. 

2. Fortuner R, Williams KJO (1975) Review of literature on Aphelmchoides besseyi 
Christie, 1942, the nematode causing "White tip" disease in rice. Helmintho- 
logical abstracts Series B, Plant Nematology 44: 1^0. 

3. Rajan AL, Mathur VK, Lai A (1989) Improved culturing technique for 
Aphelmchoidc.^ besseyi by inducing anhydrobiosis. Indian J NcmatologN' 19: 10—13. 

4. Duncan L\V, Morns Al (200G) Migratory cndoparasitic nematodes. In: Perry 
RJ\', .Morns .M, editors. Plant Nematology: CABI publishing, pp. 123-152. 

5. da Silva (iS (1992) White tip and national rice production. Informe 
Agropecuario Belo Horizonte 16; 57—59. 

6. Lin MS, Ding XF, Wang ZM, Hou FM, Lin N (2004) Description of 
Aphelmchoides besseyi from abnormal rice with 'small grains and erect panicles' 
symptom in China. Rice Science 12; 289-294. 

7. Tiwari SP, Khare MN (2003) White tip caused by Apheknckoides besseyi^ an 
important seed borne disease of rice. In: Trivrdi PC, rditors. Advances in 
Nematology. Scirntific Publishers- pp. 103-114. 

8. Togashi K, Hoshino S (2001) Distribution pattern and mortality of the white tip 
nematode, Aphelenchoides besseyi (Nematoda; Aphelenchoididae), among rice seeds. 
Nematology ?yAl-2^. 

9. Wang KM, Tsay I T, Lin YY (1993) The occurrence o{ Aphelmchoities besseyi on 
strawberry and its ecology in Taiwan. Plant Prot Bull 35: 14—29. 

10. Yu PC, Tsay TT (2004) Occurrence of a foliar nematode disease of fern in 
Taiwan. Plant Pathol Bull 13; 35-44. 

11. Tsay TT, Cheng YH, Chen HY, Lin YY, Wu WS (1995) Occurrence and 
control of nematode diseases on bulbous flowers. Plant Pathol Bull 4; 180-192. 

12. Van Megen H, Van den Elsen S, Holterman M, Karssen G, Mooyman P, et al. 
(2009) A phylogenctic tree of nematodes based on about 1200 full-length small 
subunit ribosomal DNA sequences. Nematology 11: 927-950. 

13. Kikuchi T, Cotton JA, Dalzell JJ, Hascgawa K, Kanzaki N, ct al. (2011) 
Genomic insights into the origin of parasitism in the emerging plant pathogen 
{Bursaphdenchus xylophilus). PLoS Pathog 7: el002219. 

14. Haegeman A, Joseph S, Gheysen G (201 1) Analysis of the transcriptome of the 
root lesion nematode {Pratylenchus cojfeae) generated by 454 sequencing 
technology. Molecular & Biochemical Parasitology 178; 7-14. 

15. Cheng X, Xiang Y, Xie H, Xu CL, Xie TF, et al. (2013) Molecular 
Characterization and Functions of Fatty Acid and Retinoid Binding Protein 
Gene (Ab-far-l) in Aphelmchoides besseyi. PLoS ONE 8: e6601 1. 

16. Kikuchi T, Cock PJA, Helder J, Jones JT (2013) Characterisation of tiie 
transcriptome n{ Aphelenchoides hesseyi and identification of a GHF 45 cellulase. 
Nematology 16: 99-107. 

17. Jones JT, Haegeman A, Danchin EGJ, Gaur HS, Helder J, et al. (2013) Top 10 
plant-parasitic nematodes in molecular plant pathology. Molecular Plant 
Pathology 14: 946-961. 

18. Schuster SC (2008) Next-generation sequencing transforms today's biology. Nat 
Mediods 5; 16-18. 

19. Cantacessi C, Gasser RB, Stiube G, Schnieder T, Jex AR, et al. (201 1) Deep 
insights into Dictyocmlus viviparus transcriptomes provides unique prospects for 
new drug targets and disease intervention. Biotechnol Adv 29; 261—271. 

20. Cantacessi C, Jex AR, Hall RS, Young ND, GampbeU BE, et al. (2010) A 
practical, bioinformatic workflow system for large data sets generated by next 
generation sequencing. Nucleic Acids Res 38: el71. 

21. Cantacessi C, Mitreva M, CampbeU BE, Hall RS, Young ND, et al. (2010) First 
transcriptomic analysis of the economically important parasitic nematode, 



nents. Table SI 7, RNAi inhibitors. Table SI 8, Nuclear effectors. 
Table SI 9, Neuropeptide genes encoding insulin-like peptides. 
Table S20, Neuropeptide genes encoding FMRFamide-related 
peptides. Table S21, Neuropeptide genes encoding non-insulin, 
non-FLP peptides. Table S22, Peptidases identified in A. besseyi. 
Text SI, Detailed methods of RNAseq libraries construc:tion. Text 
S2, MIRA assembly yielded sequences were not found in our 
dataset. 
(DOC) 

Files S2 Supplementary Dataset. The comparison of 

unigenes with another A. besseyi transcriptome. 

PCLSX) 

Autlior Contributions 

Conceived and designed the experiments: FW DL ZW. Performed the 
experiments: DL FW AD. Analyzed the data: FW DL. Contributed 
rcagents/matcrials/analysis tools: FW DL LL BW QC XL. Wrote the 
paper: FW DL QC. 



Trichostrongylus colubriformis, using a ncxt-gcncration sequencing approach. Infect 
Genet Evol 10: 1199-1207. 

22. Cantacessi C, Mitreva M, Jex AR, Young ND, Campbell BE, et al. (2010) 
Massively parallel sequencing and analysis of the Necator americanus transcrip- 
tome. PLoS Negl Trop Dis 4: e684. 

23. Choi YJ, Ghcdin E, Berriman M, McQuillan J, Holroyd N, et al. (201 1) A deep 
sequencing approach to comparatively analyze the transcriptome of lifecycle 
stages of the filarial worm, Brugia malayi. PLoS Negl Trop Dis 5: el409. 

24. Nicol P, Gill R, Posu-Nyarko J, Jones MGK (2012) de novo analysis and 
functional classification of the transcriptome of the root lesion nematode, 
Pratylenckus thomei, after 454 GS FLX sequencing. International Journal for 
Parasitology 42: 225-237. 

25. Yoshida K, Hasegawa K, Mochiji N, Miwa J (2009) Early embryogenesis and 
anterior-posterior axis formation in the white-tip nematode Aphelmchoides bessepi 
(Nematoda: Aphelenchoididae). Journal of Nematology 41: 17—22. 

26. Wang F, Wang ZhY, Li DL, Chen QL (2012) Identification and characteriza- 
tion of a Bur.saphelenchus xvlophilm (Aphelenchida: Aphelenchoididae) thermo- 
toierance-related gene: Bx-IISP90. IntJ Mol Sci 13: 8819-8833. 

27. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, et al. (201 1) FuU- 
length Transcriptome assembly from RNA-Seq data without a reference 
genome. Nat Biotechnol 29: 644-652. 

28. Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B (2008) Mapping and 
quantifying mammalian transcriptomes by RNA-Seq. Nat Methods 5: 621-628. 

29. Martin J, Abubucker S, Heizcr E, Taylor CM, Mitreva M (2012) Nematode.net 
update 2011: addition of data sets and tools featuring next-generation 
sequencing data. Nucleic Acids Res 40: D720-D728. 

30. Kanehisa M, Goto S (2000) KEGG: Kyoto Encyclopedia of Genes and 
Genomes. Nucleic Acids Res 28: 27-30. 

31. Moriya Y, Itoh M, Okuda S, Yoshizawa AC, Kachisa M (2007) KAAS: an 
automatic genome annotation and pathway reconstruction server. Nucleic Acids 
Res 35 (suppl 2):W182-W185. 

32. Park BH, Karpinets TV, Syed MH, Leuze MR, Uberbacher EC (2010) 
CAZymes Analysis Toolkit (CAT): Web-service for searching and analyzing 
carbohydrate-active enzymes in a newly sequenced organism using CAZy 
database. Glycobiolo^g\' 20: 1574-1584. 

33. Haegeman A, Mantelin S, Jones JT, (iheysen Q (2012) Functional roles of 
efiectors of plant-parasitic nematodes. Gene 492: 19-31. 

34. Haegeman A, Banters L, Kyndt T, Rahman MM, (ihcysen G (2013) 
Identitication of candidate effector genes in the transcriptome of the rice root 
knot nematode Meloidogyne graminicola. Molecular Plant Pathology 14: 379-390. 

35. Petersen TN, Brunak S, Von Heijne G, Nielsen H (2011) SignalP 4.0: 
discriminating signal peptides from transmembrane regions. Nature Methods 8: 
785-786. 

36. Krogh A, Larsson B, Von Heijne G, Sonnhammer ELL (2001) Predicting 
transmembrane protein topology w'ith a hidden markov model: apphcation to 
complete genomes. J Mol Biol 305: 567-580. 

37. Dalzell JJ, McVeigh P, Warnock ND, Mitreva M, Bird DM, et al. (201 1) RNAi 
Effector Diversity in Nematodes. PLoS Negl Trop Dis 5: ell76. 

38. Rawhngs ND, Morton FR (2008) The MEROPS batch BLAST: a tool to detect 
peptidases and their non-peptidase homologues in a genome. Biochimie 90: 
243-259. 

39. De Boer JM, Yan Y, Smant G, Davis EL, Baum TJ (1998) In situ hybridization 

to messenger RNA in Het£rodera glycine.s. ^ Nematol 30: 309—312. 

40. Bakhetia il, Urwin PE, Atkinson HJ (2007) QPCR analysis and RNAi define 
pharyngeal gland cell-expressed genes of (Heterodera glycines^ required for initial 
interactions with the host. Mol Plant-Microbe Interact 20: 306^312. 



PLOS ONE I www.plosone.org 



9 



March 2014 | Volume 9 | Issue 3 | e91591 



White Tip Nematode Transcriptome 



41. Kikuchi T, Jones JT, Aikawa T, Kosaka H, Ogura N (2004) A family of glycosyl 
hydrolase family 45 ccUulases from the pine wood nematode Bursaphelenchus 
xyhphiks. FEES Letters 572: 201-205. 

42. Davis EL, Hussey RS, Mitchum MG, Baum IJ (2008) Parasitism proteins in 
nematode-plant interactions. Current Opinion in Plant Biology 11: 360-366. 

43. Elling AA, Mitreva M, Gai XW, Martin J, Recknor J, et al. (2009) Sequence 
mining and transcript profiling to explore cyst nematode parasitism. BMC 
Genomics 10:58. 

44. Roze E, Hansc B, Mitre\ a \1, X'anholmr B, Bakker J, et al. (2008) Mining the 
sccrctome of the rool-knoi nematode Meloidogyne chitwoodi for candidate 
parasitism genes. Molecular Plant 9: 1—10. 

45. Jacob J, Mitreva M, Vanholme B, Gheysen G (2008) Exploring the 
transcriptome of the burrowing nematode Radopholus similis. Mol Genet 
Genomics 280: 1-17. 

46. Rosso MN, Jones JT, Abad P (2009) RNAi and functional genomics in plant 
parasitic nematodes. Annu Re\' Phytopathol 47:207-2!'^2. 

47. Tylgal T. X'anholmc B, De McintorJ, Claeys M, Couvreur M, et al. (2004) A 
new class of ubiquitin extension proteins secreted by the dorsal pharyngeal gland 
in plant parasitic cyst nematodes. Mol Plant-Microbe Interact 17: 846—852. 

48. Rosso MN, Dubrana MP, CimboHni N, Jaubert S, Abad P (2005). Application of 
RNA interference to root-knot nematode genes encoding esophageal gland 
proteins. Mol Plant-Microbe Interact 18: 615-620. 

49. Geldhof P, Visser A, Clark D, Saunders G, Britton C, et al. (2007) RNA 
interference in parasitic helminths: current situation, potential pitfalls and future 
prospects. Parasitolo^gy 134: 609-619. 

50. Bakhetia M, Charlton \VL, Urwm PE, McPherson MJ, Atkinson HJ (2005) RNA 
interference and plant parasitic nematodes. TRENDS in Plant Science 10: 362— 
367. 

51. Li C, Nelson LS, Kim K, Nathoo A, Hart AC (1999) Neuropeptide gene families 
in the nematode Caenorhabditis ekgans. Ann NY Acad Sci 897: 239-252. 

52. Li C, Kim K (2008) WormBook. Sep 25:1-36. doi: 10.1895/wormbook.l. 142.1 

53. Henrissat B (1991) A classification of glycosyl hydrolaises based on aunino-acid- 
sequence similarities. BiochemJ 280: 309—316. 



54. Lynd LR, VVeimer PJ, Van Zyl WH, Prclorius IS (2002) Microbial cellulose 
utilization: fundamentals and biotechnology. Microbiol Mol Biol Rev 66: 506— 
577. 

55. Fu Zh, Agudelo P, Wells CE (2012) Differential expression of a p-1,4- 
Endoglucanase induced by diet change in the foliar nematode Aphelmchoides 
Jragarim. Phytopathology 102:804-811. 

56. KiJcuchi T, Shibuya H, Jones JT (2005) Molecular and biochemical 
characterization of an endo-beta-l,3-glucanase from the pinewood nematode 
Bursaphelenchus xylophilus acquired by horizontal gene transfer from bacteria. 
BiochemJ 389:117-125. 

57. Karim N, JonesJT, Okada H, Kikuchi T (2009) Analysis of expressed sequence 
tags and identification of genes encoding eell-wall-degrading enzymes fi"0m the 
fungivorous nematode Aphelenchus avenae. BMC Genomics 10:525. 

58. JonesJT, Furlanetto C, Kikuchi T (2005) Horizontal gene transfer from bacteria 
and fungi as a driving force in the evolution of plant parasitism in nematodes. 
Nematology 7:641-646. 

59. Lin B, Zhuo K, Wu P, Cui R, Zhang LH, et al. (20 1 3). A novel effector protein, 
MJ-NULGla, targeted to giant cell nuclei plays a role in Meloidogyne javanica 
parasitism. Mol Plant Microbe Interact 26: 55-66. 

60. Abad P, Gouzy J, Aury JM, Castagnone-Sereno P, Danchin EGJ, et al. (2008) 
Genome sequence of the metazoan plant-parasitic nematode Meloidogyne incognita. 
Nature biotechnology 26: 909-915. 

61. De MeutterJ, Tytgatt T, Witters E, (iheysen G, Van Onckelen H, et al. (2003) 
Identification of cytokinins produced by the plant-parasitic nematodes Ileterodem 
schachtii and Meloidogyne incognita. Molecular Plant Pathology 4: 211-211 . 

62. Itiial N, Recknor J, Nettleton D, Hearnc L, Maicr T, et al. (2007). Parallel 
genome-wide expression profiling of host and pathogen during soybean cyst 
nematode infection of soybean. Molecular Plant-Microbe Interactions 20: 293- 
305. 

63. Hewezi T, Howe PJ, Maier TR, Hussey RS, Mitchum MG, et al. (2010) 
Arabidopsis spermidine synthase is targeted by an effector protein of the cyst 
nematode Heterodera schachtii. Plant Physiol 152: 968—984. 



PLOS ONE I www.plosone.org 



10 



March 2014 | Volume 9 | Issue 3 | e91591 



